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ABSTRACT 


This report traces square root information estimation, 
starting from its beginnings in least-squares parameter esti- 
mation, Special attention is devoted to discussions of the 
sensitivity and perturbation matrices, computed solutions and 
their formal statistics, consider-parameter s and consider- 
covariances, and the effects of a priori statistics. The 
constant-parameter model is extended to include time-varying 
parameters and process noise, and the error analysis capa- 
bilities are generalized. Efficient and elegant smoothing 
results are obtained as easy consequences of the filter 
formulation. 

The value of the techniques discussed here is demon- 
strated by the navigation results that were obtained for the 
Mariner Venus-Mercury (Mariner 10) multiple-planetary space 
probe and for the Viking Mars space mission. 
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SEQUENTIAL LEAST-SQUARES 
USING ORTHOGONAL TRANSFORMATIONS 


I. INTRODUCTION 

Outer space navigation analyses and o?.bit determination has relied, in the 
main, upon the method o£ least- squares introduced by Gauss (Refs, 8 and 12). 

In recent years science experiments and mission requirements have become 
more stringent, and space scientists have found it necessary to introduce more 
precise and sophisticated models, including stochastic process noise effects. 
Software specialists 5 '' were at first nonplussed because the new models 
involved time-varying parameters and process noise in a way which appeared 
to be at odds with the classical constant parameter estimation software then 
in use. 

There was a reluctance to abandon the least-squares methods in favor of 
the Kalman filter, even though the latter technique is flexible enough to accom- 
modate time-varying models with process noise. Reasons for this reluctance 
included cost, reliability and inertia. The constant parameter estimation soft- 
ware that was already developed, checked out, and proven is only a part of the 
lengthy and complex orbit determination (OD) process. In this framework the 
constant parameters are referenced to an epoch time, and thus introducing a 
current state Kalman filter would require costly modification of the entire OD 
process. Furthermore, OD problems generally involve processing thousands 
of data points and, for such situations, the Kalman filter is costly to operate. 
Least-squares analyses using or* *ogonal transformation techniques (Ref. 10) 
have proven to give reliable and accurate results, while the Kalman filter in 
process noise free situations exhibits numeric deterioration and instability. In 
light of these comments, it is not surprising that software specialists and man- 
agers were reluctant to abandon least-squares. 

The apparent contrariety between the least- squares parameter estimation 
techniques based on orthogonal transformation methods and the Kalman covari- 


The remarks and opinions which follow reflect the author's experiences at the 
Jet Propulsion Laboratory but are believed to have a more general applicability. 
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ance oriented methods was cleverly resolved by modifying both the model and the 
least-squares orthogonal transformation method, cf Curkenda.il and Leondes 
(Ref. 5) and Dyer and McReynolds (Ref. 6). The Dyer -McReynolds filter algo- 
rithm, dubbed with the acronym SRIF (Square Root Information Filter) by 
Kaminski (Ref. 9 ) » is of primary importance in our estimation developments 
and discussions. The work in this report documents estimation techniques used 
at the Jet Propulsion Laboratory which have been found to be efficient, have 
good numeric characteristics, and which have led to new and improved insights 
into the problem of orbit determination. 

The intention here is to summarize and highlight the epoch-state model 
(see Eq. 3-4) SRIF formulations for filtering, smoothing, and error analysis. 
Specializing the results to the epoch-state model greatly simplifies the some- 
what intimidating algorithms of Refs, 2, 3, 4, and 9 that were developed for 
more general linear models. Our techniques and algorithms arc applicable to s, 
rather broad class of problems even though our development is oriented toward, 
and was motivated by, orbit determination applications, The following paragraphs 
outline the scope of material that is to be discussed, 

In Section II we discuss Golub's sequential solution to the least-squares 
parameter estimation problem through the use of Householder orthogonal trans- 
formations* (Ref. 7). This method is acknowledged to be numerically stable and 
is considerably more accurate than the classical normal equation technique intro- 
duced by Gauss (Ref. 10), Various jargon are defined here: the "data equation" 
corresponding to an estimate -estimate error covariance pair; "computed" solu- 
tions and covariances corresponding to restricted dimension models; "sensitivi- 
ties" and "perturbations" corresponding to neglected parameters, and "consider" 
covariances which reflect the true error covariance of the computed estimate. 

In Section III the "epoch state" model is defined along with a technique for 
including time propagation and white process noise. Certain computational 
efficiencies and algorithmic simplifications come about when the process noise 
entering the time propagation is exponentially correlated. 

Smoothing, discussed in Section IV, is arranged so that it blends with the 
material of Sections II and III. Smoothed analogues are developed for the 


Householder orthogonal transformations are discussed in Appendix A. 
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"computed" estimate, the "computed" estimate error covariance, and the 
"sensitivity" matrix, These results differ from their filter counterparts, which 
arise from auxiliary information array computations, in that they are computed 
recursively. Commonalities between the filter and smoother algorithm formu- 
lations allow for a generalization of the "consider" covariance and of estimate- 
covariance pairs corresponding to various order estimation models. 

Section V is devoted to error analyses that are related to the SRIF. The 
data equation representations lead straight-away to augmented information 
array algorithms. The error equation methodology parallels the covariance 
error analysis methods that are derived in Ref, 11; but the error analysis algo- 
rithms, like the SRIF, differ from their conventional covariance algorithm 
counterparts. The error analyses are not completely general but do cover the 
principal error sources occurring in orbit determination, including the effects 
of incorrect a priori state and process noise covariance statistics, unestimated 
parameters, and the effects of incorrect colored noise time constants. 


preceding page 
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II. 


SEQUENTIAL DATA PROCESSING USING ORTHOGONAL 
TRANSF ORMATI ONS 


Accurate and stabLe sequential least-squares data processing methods arc 
of major importance in orbit determination because OD problems generally have 
poor observability characteristics (are ill-conditioned) and involve large amounts 
of data. The normal equation methods, classically used in such problems, have 
been supplanted by the orthogonal tx*ansformation techniques that were introduced 
by Golub (Ref, 7), To see how orthogonal transformations enter into the least- 
squares problem, consider the overdetermined >, ‘ system: 


z = A x + v (2-1) 

where x is an n-vector, z is an m-vector (with m > n) and v is assumed to be 
scaled so that vcN(0, I m ), The least-squares problem is formulated and solved 
in the following sequence of equal-ons; 



Ax - z 


2 


= ||Q a (Ax-z )|| 2 


( 2 - 2 ) 

(2-3) 

(2-4) 


Q. such that Q A Q A = I 
A A A m 


Q a Ax- Q a z 


q a a 


(2-5) 


|m-n 


"A 

Z 


, Q a z = 




m-n 


' Only overdetermined problems are discussed here because when a priori infor- 
mation is used the estimation problem is of this form. 
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j £s ( x ) « 


(2-6) 





(2-7) 


Golub pointed out that the orthogonal transformation could be simply com- 
puted using Householder's reflections, that in this ci.ie R will be a triangular 
matrix, and that this reduction of the array [A, z| could be accomplished without 
explicitly computing Qy^. Appendix A reviews the key ideas and formulae 
involved with Householder matrix reductions I v 

From Eqs. 2-2 and 2-7 it is evident that J^ g is minimized for x chosen 
such that 



( 2 - 8 ) 


and that ||e|| is the minimum sum of squares residual error. Statistical signifi- 
cance for the solution x is an easy consequence of Eq, 2-1, i, e. , 


Q a z = C 

I a A x + Q a 

V 

A _ 

1 Z 


_ A 

R 


_ A 
V 

I e ! 


L° . 

X + 

V 

L e 


and thus, 

A A A 
Z s Rx + V 


(2-9) 


where vtN(0, I n ), because N (°» Combining EqB. 2-8 and 2-9 gives 


A 

X - X 



I 


‘For more details of this subject the reader is referred to the book by Lawson 
and Hanson (Ref. 10) which comprehensively treats the various numeric and 
computational aspects of orthogonally computed least-squares solutions. 
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so that 


Pa a E [(x - x) (x - x) T ] = R" 1 R“ T (2-10)* 

A _ 1 

Equation 2-10 shows that R is a square root of the error covariance matrix PA, 

A 

An equation of the form of Eq. 2-9, with R nonsingular, is called a data 
equation . From Eqs. 2-8 and 2-10 it is clear that there is an equivalence 
between the Information array [A, z] and the covariance -estimate pair (Pa, x). 
Information arrays are not unique, but if [Rj, Zj] and [R^, z,,) correspond to the 
same covariance -estimate pair, then for some orthogonal matrix Q, 

[*l,*l3 = Q Z z^ * 


(This relation Is a direct consequence of the data equation representation corre- 
sponding to the information array. ) An obvious but important consequence of a 
data equation representation for the a priori covariance and estimate is that the 
least- squareB orthogonalization is recursive. This result is written as 
theorem 1, 

Theorem 1 (Recursive Least-Squares Measurement Updating) 

T 

The data z. = A.x + y., with E(v.) = 0, E(y., v . ) = can be recursively 

JJJ J rAA^ 

least-squares fit as follows. Start with the a priori information array LRq, Zq], 

and for j = 1, . • * , N, recursively compute: 


PA 

A ~ 


A 

A “ 

R j-> 

z j-i 


R. 

Z. 


J 

J 

A. 

Z. 


0 

e. 

J 

J 

J 




2 = 
j 


2 + 
j-1 + 



• where 



orthogonal, is implicitly defined by Eq. 2-11. 



( 2 - 11 ) 


( 2 - 12 ) 
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Then, for each j, Xj given by Eq. 2-13 is the least-squares estimate of x based 
on the data {zji 55 0, » • • , j}, and Pj is its error covariances 


A 

X. 


A-l A 

= R j "j 


A A - 1 A _T 

P, = R, R. 


(8-H) 


Moreover, 


2 




(2-14) 


NOTES: 

A 

(1) R. are by construction triangular matrices, generally arranged so 
J 

as to be upper triangular. The triangular structure reduces compu- 
tations in Eq. 2-11 and also facilitates the matrix inversion appear- 
ing in Eq. 2-13. 

(?) The recursion (Eq. 2-11) is similar in certain respects to the 
normal equation results, 

[vj ■ [ A j-i'V>] + [ A 7 V\ T z j] < 2 - 15 ' 

In this case the estimate and covariance are given by 

*j = A j‘ ! d j ■ * A j' (2 - 16) 

Equation 2-15 is susceptible to roundoff errors; and not infrequently 
one finds that the computed A is singular, cf Refs. 7 and 10. 

(3) Estimates and covariances (Eq. 2-13) need not be computed for each 
j and this allows for certain economies of computation. 

(<) Equation 2-12 gives the accumulated residual error corresponding 
to the latest estimate; but this computation does not explicitly 
involve the estimate. 
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(5) The relations in Eqs. fc-i c ar.d 7 are important in residual 

analyses problems involving large amounts of data. In such prob- 
lems it is inconvenient and costly to reread the previous data 
{A., zji = 1, • • • , j [from storage and compute Eq, 2-14, 

The data equation, as shown in Theorem 1, is a handy vehicle with which 
to express the state of the system. To further explore this representation, 
suppose that x is replaced by the partitioned vector and that the data Eq. 2-9 
is partitioned consistently. Assuming that R is upper triangular gives: 


— 

— 


— — 


— — . 


— — ■ 

R 

R 


X 


z 


V 

X 

xy 




X 


X 

0 

R 


y 


z 


V . 


y 




y 


y 


(2-17) 


x and y are now assumed to be vectors of dimension n x and n^. respectively, 

/ith n = n + n . 

x y 

interpretation. 


This partitioned data equation ha,s an interesting and useful 


The bottom equation, R^y = z^ - v is a data equation for y and as such it 
implicitly contains the estimate and covariance of y. This point is of importance 
when the parameter estimation problem is of large dimension: but one is oi.^y 
interested in explicitly estimating a small subset of parameters. 


Turning attention to t-he upper of Eqs. 2-17, we point out that a 
consequence of the triangular construction is that R and v are independent 
of y. This observation is of major importance because it enables one to 
analyze the effects of various y models in a straightforward and insightful 
manner. With these comments in mind, Eq. 2-17 is expanded to obtain: 


x = R" 1 z + ( -R" 1 R ) y - R" 1 
xxyx xy / 7 x 


= x c + (SEN) y - R x v x 


(2-18) 


This innocent looking equation is actually pregnant with implications, which 
we now proceed to describe. 
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X 


c 



z 


X 


SEN = -R - 1 
x 




( 2 - 19 ) 


where x is called the computed estimate and corresponds to a model with no 
y parameters ; SEN is called the Sensitivity ; and called the computed covariance 
will equal E[(x - x c ) (x - x c )T] if the model contains no y parameters. 



{ 2 - 20 ) 


where P is the a priori covariance of the y parameter vector. P Xc is the 
actual or true covariance of the estimate x c ; it reflects the estimate error due 


to using a filter model that ignores the y parameters. 


Not infrequently one includes y parameters in the model, not because they 
are of special interest, but because excluding them may lead to poor estimates 
of the x parameters. Equation 2-20 separates out that portion of the error 
covariance that is due to the omitted y parameters. In this context the compo- 
nents of y are called consider parameters and P Xc is called a consider 
covariance. 


In the special but important case where the y a priori covariance is diag- 
onal, one can get further useful information about the effects of the individual 
components of y. In this case, the columns of JT s SEN Diag (oy(l) f • .o^ny)) 
are, so to speak, error perturbations of the estimate error; and perusal of qis 
helpful in determining which y parameters contribute most or least to uncertain- 
ties in particular components of the estimate. This information is a valuable 
aid in identifying major error sources. 


The importance and utility of the SEN matrix does not end here; it appears 
in the optimal least-squares estimate and covariance results. From Eq. 2-18 
it follows that 


x = x +SENy 
c 


(2-21) 
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E [(x - x) (x - x) T ] 


( 2 - 22 ) 


where y = R 


-1 

y 


z and 

y 



= P__ + SENP y (SEN) T * 


x 


a -1 -T 

where P = R R 

y y y 


Furthermore, from Eqs, 2-18 and 2-21 it follows that 


SEN 


9x _ 9(x - x) 
3y 3(y - y) 


(2-23) 


and thus SEN is a sensitivity matrix which measures the sensitivity of the esti- 
mate to the y parameters. 

These formulae involving the sensitivity matrix point up its key roLe in 
parameter estimation analysis. Although these estimation formulae are easily 
and simply implemented, their flexibility and utility is enhanced by noting the 
following: 

(1) Systematically increasing the x parameter set can be done with but 
a modicum of effort. The necessary equations follow from the 
accompanying matrix identity: 




-i 


(-R~ l R ^ R ’ 1 

R 

R 


X 

xy 


X 

\ x xy/ y 

0 

R 


0 

R" 1 


y J 



y J 


( 2 ) 

( 3 ) 

( 4 ) 


SEN is independent of the y a priori and this simplifies certain 
analyses. 


The triangular structure of R^ makes it a simple matter to delete 
the lowermost y's and examine models with a smaller number of 
parameters. 


One need not commit to a prespecified arrangement of the y param- 
eters; i. e. , it is possible to omit any subset of parameters from the 
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model. In the general case, one should permute the columns of R^ 
and SEN corresponding to the desired y-parameter arrangement 
and then use an orthogonal transformation to retriangularize R^. 


To illustrate the kind of analysis that can be performed using the ideas 
just described, consider the following hypothetical but typical problem. Formu- 
lation of the problem involves, say, 150 parameters; but we are interested in 
estimates for perhaps 10 of these, and they are designated x. Most of the 
remaining parameters came about because of our efforts to precisely model the 
problem. Based on our knowledge of this problem, we arrange the parameter 
name list so that those parameters believed to be most significant are at the top 
of the list. There are 3000 data points and these are orthogonally transformed 
to an upper triangular form, cf Eq. 2-11, The sensitivity matrix is formed, 
cf Eq. 2-18, and is 10 by 140, A perturbation matrix, F, is computed using 
nominal y-parameter a priori standard deviations, and the columns are 
inspected. This suggests a rearrangement of certain of the y's, i. e. , moving 
to the end those y's that appear to have negligible significance. The columns of 
SEN and Ry. are rearranged accordingly, and Ry is orthogonally triangularized 
and the y a priori is included. It is possible now to exhibit estimates and esti- 
mate error covariances based on, say, 50, 60 and 70 y-parameter models. If 
the higher order y's are truly less significant, the estimate and covariance 
should be relatively unaffected by increasing the model order. By redefining 
the x and y namelists, or by partitioning and expanding the appropriate equations, 


one could examine the effects of estimating the first n^ parameters, considering 
the effects of the next ny-, parameters and ignoring the rest. For a problem 
with this structure, we would be perusing P x _ where xj and y^ represent 

rp C C 

(* , Yp * ’ '* Y n ) and (y n +1 , • • • , y +n ) respectively, 

y i Yi y i y z 

The data analysis options described in this section are standard compo- 


nents of JPL's Orbit Determination Program. 
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III. MAPPING AND THE INCLUSION OF PROCESS NOISE 


The discussions of the previous section were directed at constant 
parameter estimation. Fortunately, and not coincidentally, most of the key 
results and notions of that section extend to the time varying' problem, 



4>.X. + G.w. 
J J J J 


(3-1) 


where ($.} are nonsingular and {w.} is a white noise sequence. The nature of 
J J 

our orbit determination problems motivates us to partition X, and we write 


p 


‘ M 

0 

0 ' 


p’ 


w j 

x 

= 

d> 

T xp 

4> 

r x 

Ky 


X 

+ 

0 

.y. 

j+i 

0 

0 

I- . 

j 

. y. 

j 

_ 0 _ 


where the 4> elements are transition matrix elements <b (j+l,J), etc. The 
vector p is colored process noise and y are constant parameters. In this 
formulation M, assumed to be diagonal, may be singular; and thus, we allow 
for white noise inputs to x. The model of Eq. 3-2, while somewhat special- 
ized, is quite adequate for our applications. Spacecraft related p parameters 
are generally random accelerations; viz solar pressure, altitude control jet 
gas leaks, or satellite drag effects. OD related p parameters may reflect 
ephemeris variations or Deep Space Network Station location position 
uncertainties. 

Previous orbit determination analyses which did not involve process 
noise were cast as parameter estimation problems by referencing the time- 
varying parameters to an epoch time. We retain this practice because such a 
formulation is compatible with existing OD software (at JPL) and because our 
estimation computations are somewhat simplified in this framework. Thus, we 
define the epoch state variable x., so that 


Xj = 4> x {j* °) x j + 4> xy (j* °)y 


(3-3) 


x. are called epoch state variables because when 4>..„(j + i, j)p ; — 0, x. reduces 

J _ X P J J 

to the epoch state Xq. When Xj is substituted into Eq. 3-2, one obtains the 

simplified propagation representation (Eq. 3-4), 
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’ p" 


’ M 

0 

0‘ 


w m 

P 

* 

* w." 

J 

X 

= 


I 

0 


X 

+ 

0 

„y. 

j+ 1 

0 

m 

0 


j 

.y. 

j 

. o. 


where V {j) = 4>“ l {j+l, °) 4> xp (j+ l * J)* 

Let us turn now to the business of constructing the SRIF mapping algorithm 

corresponding to the dynamic model (Eq. 3-4). The mapping process consists 

simply of constructing a data equation for given an a priori data equation 

for X.. Thus, suppose we are given an a priori data equation in the following 
J 

form; 


A 

R 

P 

r\ 

R 

px 

r 

o. 


' p j' 


_ A „ 
z 

P 


_ A 

P 

A 

R 

xp 

A 

R 

X 

A 

R 

xy 


X. 

J 

= 

A 

z 

X 

- 

A 

V 

X 

o 

1 

0 

A 

R 

y J 


.y. 


A 

z 

L yJ 


A 

L "y J 


NOTES: 


(3-5) 


(1) In the interest of notational economy we omit the j subscript on R 
and z, and rely on the symbol "a" to indicate the values of R and z 
at the start of the mapping. 


(2) A pragmatically important advantage accruing to a formulation 

with the p variables uppermost is that when R = 0, the formulae 

xp 

of Section II are valid and need no modifications. 


Using Eq. 3-4 it is oasy to replace x^ by Xj +1 , i* e. 


- 



■ 


- ■ 


- • 


- 

A 

R - 
P 

A 

R V 
px p 

A 

R 

px 

A 

R 

py 


p j 


A 

Z 

P 


A 

v p 

A 

R 

xp 

A 

- R V 

X p 

A 

R 

X 

A 

R 

xy 


x j+i 

= 

A 

z 

X 

- 

A 

V 

X 

. 

0 

0 

A 

R 

y_ 


y 


A 

z 

L yJ 


A 

v y 

m 


(3-6) 


* 


Generally the R terms are the resultofdata processing, 


and in this care R = 0, 
xp 
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To complete our construction of the data equation at (j+1), wo have to 
roplacc pj by Pj + i» and this is done in the following way. Recall from Eq, 3-4 
that the colored process noise terms satisfy 


p . , , = Mp, + W„ 

l j+l ‘j / 


(3-7) 


and that the statistics of w. can be written in data equation form 

J 


R W . - 7 , - \> 

W J W W 


(3-8) 


(z =0 wb<sn E(w.) = 0. ) 

' w ' j' 

Equations 3-6 through 3-8 are combined by considering an augmented 

state vector, one that contains p. and p. 

J J 


f - R M 
w 

R 

w 

0 

0 ‘ 


— m 

P J 


1 

£ 

<53 

1 


N 1 

w 

A A 

r - n v 

p px p 

0 
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R 

px 

A 

R 

py 


P j+1 


A 

z 

p 


A 

V P 

A A 

R - R V 
xp X p 

0 

A 

R 

X 

A 

R 

xy 


*j+l 


A 

z 

X 


A 

V 

X 

0 

0 

0 

A 

R 

y J 


y 


1 

N> 

1 


.V 


(3-9) 


To obtain the desired data equations from this result one has only to orthogonally 

eliminate p. from the bottom equation; and, in information array form, the 
J 

result is: 
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Pj 


Pj+i 
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( 2 - 10 ) 

where Q, orthogonal, is chosen to partially triangularize the information array 
and the S notation represents the time updated information array 


n 


P 




NOTES: 



(3-11) 


(1) The top row of Eq. 3-10 corresponds to a data equation 


( 2 ) 


R p. + R p.., + R x., . + R y - z. - v , 

P*J PP\J+* px j+l py y P P 

and this equation is of major importance in smoothing, 

Section IV, 


(3-12) 

cf 


The y information array J z^J is unchanged by the mapping. 

(This is as it should be since the y parameters are time- 
invariant. ) 


(3) Appendix B contains an efficient and compact FORTRAN mech- 
anization of the mapping (Eq. 3-10), This mechanization makes 
maximal use of the special structure of our problem formulation. 

(4) Dyer and McReynolds (Ref. 6) describe how to time update an 
information array corresponding to the general dynamic 
model (Eq, 3- 1) , 
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IV, SMOOTHING 


Most orbit determination involves nonreal time or postflight data 
reduction and, as a result, there is an interest in smoothing, l.o,, trajectory 
estimation based on the entire data history. In this section we show that 
smoothed estimates can bo computed with but a modicum of computational 
cost. We also present smoothed analogues of the computed estimate, computed 
covariance, the sensitivity matrix, and the consider covariance. These gener- 
alizations are important because the formulae are consistent with the constant 
parameter results of Section II; and thus the software developed for that 
problem can be applied to this problem as well, Furthermore, OD engineers 
familiar with constant parameter estimation can better relate to the smoothing 
problem when it is couched in a familiar framework. 

To see how the smoothed state estimates come about, suppose that 

smoothed estimates of Pj + j» ^j+1 anc * V ^ ave been computed and that these are 

labeled with a Smoothed estimates for p. and x. are a direct result of the 

J J 

mapping, Eq. 3-4, and the key portion of the data equation, Eq. 3-12, Define, 
cf Eq. 3-12, 



S|{ 

and note that the smoothed estimate satisfies the data equation with the v term 
set to zero, i. e, , 


* T T * T * T * 

p. = L - L| p.,. - L x. , , - L y 

\J z p 1 J+l x j+1 y 7 


x.’ = X. ,. - V p.’ (cf Eq. 3-4) 

J J+l P \J H 


(4-2) 

(4-3) 


Equations 4-2 and 4-3 establish the desired smoothing recursions. 
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NOTES: 


(1) 

( 2 ) 


The reader is reminded that the L torms and V 

P 


are j dependent. 


The smoothing recursion, which is backward in time, is initialized 
with the filter estimates of the terminal time point. Smoothed esti- 
mate error covariance recursions could bo obtained in the usual 
fashion (Ref, 11) by differencing Eqs. 4-2 and 4-3 with Eqs. 3-4 
and 3-12, squaring (multiplying by the transpose), and applying the 
expectation operator. 


An alternative and better formulation of the smoothing solution results if 
the algorithm is arranged in a form that is analogous to the SRIF solution for 
the filtering problem with a computed estimate and error c^ariance and a 
sensitivity matrix. With this in mind, we decompose the p-x ver+or as 


p 


Pc 


SEN* 

P 


= 


+ 

$ 

X 


X 

c 


SEN 

X 

m m 


(4-4) 


where the "c" subscripted variables are independent of y. Substituting this 
decomposition into Eqs. 3-4 and 3-12 leads to 


= L Z - L p P c <j +1) - L X x c<j + l) - R p _1 V P 


x c (j) = x c<j +1) - y P P C W 


(4-5) 


SEN* (j) = -h SEN* (j+1) - L SEN* (j+i ) - L> 
p p ' x x J ' y 


SEN*(j) = SEN* (j + 1)- V p SEN* (j+1) 


(4-6) 


The smoothed state recursion for p c and x c follow from Eq. 4-5 by setting v p 
equal to zero; 
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To obtain a covariance recursion for the p and x estimate errors, we differ- 

* c c 

once Eqs, 4-5 and 4-7? 


AP C 


Ax 


c 

m m 

j 



V L 
P P 


-JL 

x 

I+V p L x 




Ax^ 


j+l 



(4-8) 


where Ap = p - p and Ax = x - x . From this, the covariance recursion is 
1 c r c * c c c c 

seen to be: 



tig _ jij 

From the way that x' , P c ' and SEN^ have been defined, it follows that the for mu- 
lae of Section II are applicable, i. c, , 



P 


X'*’ 


P c + SEN P SEN 

x X y X 




* 1 . 

X + SEN y 
c X ' 


C ! !< ! l ! A s|t T 

P c + SEN P (SEN ) A 


(cf Eq. 

2-20) 

(4-10) 

(cf Eq. 

2-21) 

(4-11) 

(cf Eq. 

2-22) 

(4-12) 
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<cf Eq, 2-23) 


(4-13) 




SEN* a 8 ^ x " 

x a(y - y * ) 


NOTE; When p is a scalar, the filtering and smoothing computations reduce 
considerably. This observation Is exploited in the FORTRAN mech- 
anizations given in Appendix B, where the algorithms are arranged 
so that the effects of a vector p arc included one component at a t'me. 
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V. 


ANALYZING THE EFFECTS OF THE A PRIORI STATISTICS 


Engineers and scientists who choose to apply linear estimation technology 
are often perplexed about what values to use for tho a priori statistics that 
appear in the estimation algorithms. One generally has only an order of 
magnitude knowledge of tho a priori uncertainties and, because of thiB, there 
may be doubts about the reliability of the results. Examples of the kinds of 
questions that should be answered in order to intelligently apply estimation 
software to meaningful problems are: 

(1) How would doubling or halving the a priori state vector variances 
effect the estimates and estimate error covariances? 

(2) Suppose that the measurement errors were actually larger {or 
smaller) than assumed. How would this affect the estimates and 
covariances ? 

(3) Process noise and correlation times aie mathematical fictions 
designed to compensate for various modeling errors and physical 
phenomena. How sensitive is the estimation algorithm to these 
to rms ? 

Answers to questions of this nature aid in filter design and create an 
atmosphere of confidence. The effects of bias parameters arc best analyzed 
in terms of the sensitivity matrix and the consider covariance, and these have 
been discussed in Section II. In this section, we concentrate on the process 
noise effects and, to keep tho notation from becoming cluttered, we omit the 
constant y parameters from this discussion (note that their omission is only 
for expositional simplicity; the mathematics is not at all complicated by their 
inclusion). 

To analyze the effects of incorrect a priori statistics in tho data equation 

estimation formulation (Eq, 2-9) one need only maintain a knowledge of the 

covariance of the v error term. To better understand what is meant by this 

statement, and to initialize the filter error algorithm, suppose that 

X = jj^J , is assumed to have an a priori estimate error covariance, 

P = r“1 
X K X 

If we now write 


R 


-T 


but that the actual error covariance is 


p x - (*x)-‘ K)- 


v x 


= R x x 


+ 


A 

= R x x 


(5-1) 
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then will not have an identity covariance. Instead, 

E ( v x v x) ~ A x A x’ A x = R x( R x) (5_2) 

A 

because the true data equation for the estimate X is 

R“X = R“X+^, E[4(4) T ] « I (5-3) 

a 

and comparing Eq. 5-1 with Eq, 5-3 gives Thus, is a 

square-root of e ( v x v x )’ Note that one can recover the actual covariance 
from a knowledge of A^ and R^, i. e. , 


p x = ( r x 1 a x) ( r x 1 a x) 


{5-4) 


The error analysis technique for the filtering algorithm is to sequentially 
update the augmented array 


[ R x ' ' A x] 

To see how A^. is updated when data is addended suppose that 

z = A X + v 

is to be included with the a priori: 



assumed 


actual 


(5-5) 


(5-6) 


(5-7) 
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then the SRIF algorithm would 


compute 


[ 


X 



from 



(5-8) 


It is easy to show that A^, the updated square root covariance, can be 
obtained from 




. ■ 


m 

i 

■ 



! o 

i 


1 

r„ r ! 

r. 

A 

X 


x ! 

Xv 

Q 

1 

1 

1 

1 

O t 

1 
1 
i 

1 

1 

;.R (r*) 1 


I 

l 

1 

* 

i 

*•* 


L 

• v\ v) J 


i 



(5-9) 


The F array is a square root of ^( v x v x)' ^ut ^ * s rec ^ an S u ^ ar ‘ For reasons 
of storage, it is preferable to compress V with a triangular form: 


X 



(5-10) 


where Q orthogonal is chosen to triangularize the T array. It is clear that 
a a a t t <p 

obtained this way is such that A^A^ = * ^Xv^Xv" 

The importance of Eq, 5-9 is that it shows that the T terms can be 
obtained as addended columns to the filter algorithm, i. e., 
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r 


r R 


A 

Q 


X ' 

1 
1 
» 

z x : 

i 

i 

i 

*x : 

i 

i 

i 

0 

1 

1 

1 

A ' 

i 

i 

R,. z i 

i 

i 

o i 

R (R*) 

v • 

V i 

i 

v\ v 1 

A 

R x 

' A 

; z x 

i 

i 

i 

• r 

! x 

i 

i 

1 

1 r 
; Xv 

i 

i 

0 

1 

t 

1 # 

i 

i 

I # 

i 

i 

! * 


<5-11) 


The time propagation portion of the SRIF algorithm is handled in the same way 
as was depicted in Eqs, 5- 10 and 5-11, i, e. , 


- 1 -T 
R 1 R x 
w w 


assumed 


T 

E(ww ) = 


v-T 


(5-12) 


actual 


[(<)" 1 ( R w)" 

and the combined time propagation algorithm is (cf Eq. 3-10) 


R* 

1 

R* 

1 

R* 

1 


1 

l 

# 

1 

1 

* ■ 

P 

1 

1 

1 

PP 

1 

1 

1 

px 

1 

1 

1 

p 

i 

1 

i 


1 

1 

1 


0 

1 

1 

1 

1 

1 

S 

P 

1 " r 

i 

i 

i 

i 

i 

rw 

S 

X 

• 

1 

1 

1 

1 

1 

s 

y 

t 

r 

i 

i 

♦ 

i 

r 

1 X 

1 

1 

1 

1 

t 

1 

r Xw . 


Q 


r -r m j 

w ; 

i 

i 

R 

w 

! 0 

1 

1 

1 

1 

» 

z 

w 

1 

I 

1 

1 

t 

0 

1 

* 

< 

) n 
i p 

i 

• 

A A { 

s - s v ! 

p ^ p • 

0 

! a 

s 

1 X 

V 

1 

1 

1 

1 

1 

A 

s 

z 

1 

1 

1 

1 

1 

A 

A x 

1 

1 

1 

! o 

) n + n 
P x 


(5-13) 
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(5- 14) 


where A 

w 



and 


Pa ! ( 

>1 = r? ! 1 

L A x ' 

J L x Xw J 



with Q w orthogonal. 

NOTE: The entries designated with a "’!*" represent square root error 

covariance terms that contribute to the actual smoothing covariance. 

One need not derive additional equations to evaluate the effects of using 
the wrong colored noise, i. e. , using colored noise terms in the filter model 
which have different time constants from those of the actual model. One 
merely includes the variables in question twice, labeled say as and p & . The 
p^ variables assume the filter values and have zero variances’!* in the actual 
model; the p variables are assumed to have zero variances’!* in the filter 
formulation. This simple state augmentation procedure avoids burdensome 
additional software, 

A jr included in the filter algorithm as just described includes the effects 
of using an incorrect initial covariance, incorrect measurement a priori, 
incorrect process noise a priori, and mismodeling of the colored noise in the 
filter model. Studying these kinds of effects promotes insight into the filtering 
process, reveals limitations of the algorithms, and most importantly it influences 
filter design. 


’^Information filter formulations can not accomodate zero variances. In 
practice, small values (based on engineering 'judgement and computer word- 
length considerations) are used; the results have been quite satisfactory. 
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APPENDIX A 


PROPERTIES OF HOUSEHOLDER ORTHOGONAL TRANSFORMATIONS 


One could use a variety of orthogonal transformations to effect the various 
matrix triangularizations that have been referred to in this paper. In this 
appendix we discuss a class of orthogonal transformations (Householder 
reflections) which have been quite adequate for our purposes. Householder 
transformations have the form Q(u) = 1-2 puu^; p - * = u^u where the components 
of u are generally chosen to null out certain components of some predetermined 
vector; this point is amplified in properties A-l and A-2, below. 

The matrix Q(u) plays no explicit role in our estimation work. However, 
the following list of properties are key to successful application of these 
transformations . 


(A-l) 


(A-2) 


If u is chosen to zero all but the firBt component, a^, 
vector "a", then one can take u. = a-, for i > 2, and u, 
+ sgn (aj) v a"^a; 

In this case Q(u)a = ( -sgn (a^) \/a T a, 0, . . , , 0) T , 


of 


= a 


1 


If u. = 0, then Q(u) acting on any vector b (i. e. Q(u)b) leaves b. 

J J 

unchanged. 


(A-3) The matrix Q(u) is an implicit quantity; Q(u)b is expressed 

directly in terms of b and u as Q(u)b = b-Xu, where X = 2p(u^b). 

(A-4) If b is orthogonal to u then Q(u)b = b. 

T 2 

(A-5) Q(u) = Q (u) and Q (u) = I; i.e. Q(u) is a symmetric and 
orthogonal matrix. 


These properties enter into matrix triangularization of a matrix A by 
applying a sequence of these elementary transformations. At the jth stage, 
one chooses a vector whose first j-1 components are zero and which zeros 
out all elements of column j (of the (j-l)st transformed matrix A^” 1 ^) that are 
below j. The jth transformed matrix is A^ = Q(u^) A^~^. Observe that 
Q(u^) has no effect on the first j-1 columns of A^ - , because they are 

orthogonal to u^'. Also, the first j-1 rows of A^ - ^ are left unchanged because 
of property A-2. Thus, Q(u^) acts only on the lower matrix partition of A^”^ 
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composed of the last m-j rows and n-j columns, (A(m,n)), and it zeroes out 
all elements but the first in the first column of this subarray. This method is 
illustrated in the FORTRAN mechanizations given in Appendix B. 

The triangularization procedure applied to 

Ax - z-v 

where v has unit covariance, is analogous to a Gaussian elimination process that 
preserves the unit covariance characteristic of v. 
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APPENDIX B 


EPOCH STATE FILTER/SMOOTHER FORTRAN MECHANIZATION 


The filter /smoother algorithms described in Sections II - IV are FORTRAN 
mechanized in this appendix. The main reason for including FORTRAN mech- 
anizations is that in this way we can best show how compact and efficient our 
algorithms are. 


Essentially all the filtering computations take place within an array 

Sin + n +m,n + n + n +1) and a vector R n (n + 3)/2 where n , 
\px p x y / y L y ' y * J p' 

n x> n and m are the respective dimensions of p, x, y and the number of 
measurements , 


NOTE: The FORTRAN descriptions to follow do not correspond to exact pro- 

gramming conventions. Greek symbols, lower case letters, square 
root and inequality symbols, etc, , were purposely used in an endeavor 
to make the code more readable. 


I. FILTER: MEASUREMENT UPDATING 


Inputs 

S(Npx +m, Npxyl); where Npx = n +n and Npxyl = Npx + n + 1, The 

p X y 

upper Npx rows of S contain the a priori array corresponding to Eq. 3-11; 
and the bottom m rows correspond to the m measurements, each with unit 
variance. 



tains the a priori 


+ 3 

y 


)/2j, a vector stored upper triangular matrix, 
information array. 


which con- 


Outputs 

The updated SRIF information array occupies the upper Npx rows of S, 
and the vector R contains the updated SRIF information array for the y 
parameters . 
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DO 40 J = l, Npx 

cr = z @z = ze'-'o 

DO iOI=J, L @L = Npx + m 

w(I) = S(I, J) @w(Npx): work vector 

S(I, J) = z 

10 cr = or + w(I)**2 
If (or < z) go to 40 

c If cr equals zero, col, J Is zero and this 
c step of the reduction is omitted 
cr = Vff 

If [w(J) > zj cr = -cr 
S(J, J) = cr 
w(J) = w{J)-cr 

cr = one/ |V* w(J)J @one = 1. 

JP1 = J + 1 

DO 30 K = JPl, Ntot @Ntot = Npx + n + l 

6 = z 

DO 20 I=J, L 

20 6 = 6 +S(I, K)*w(I) 

6 = 6*cr 
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DO 30 I=J, L 

30 S{I, K) = S(I, K) + 6*w(I) 

40 Continue 
c 

c This completes the S portion of measurement processing 
c 

e Begin y-reLated portion of measurement processing 
c 

KK = 1 

DO 90 J = 1, Ny @Ny = n y 

a- = Ry(KK)>:c>:c2 @Ry = R^ 

DO 50 I = 1, m 

50 <r = a; + S(Npx + I, Npx + J)>!<*2 
If ( a* = z) Go to 90 
cr = \[<t 

If [Ry(KK) > z] cr = -<r 
6 = Ry(KK)- cr 
Ry(KK) = cr 
P = one/{cr*6) 

JJ = KK 
L •- J 
JP1 = J + l 
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DO 80 K = JP1, Ny l @Nyl = Ny + 1 

JJ = JJ + L 

L = L + 1 

c r = 6>!<Ry(JJ) 

DO 60 I = 1, m 

60 <r = <r + S(Npx + I, Npx + J)*S(Npx + I, Npx + K) 

If (tr = z) Go to 80 
cr s o->!<Beta 

Ry(JJ) = Ry(JJ) + <r#6 
DO 70 I = 1, m 

70 S(Npx + I, Npx + K) = S(Npx + I, Npx + K) + cr>j<S(Npx + I, Npx + J) 

80 Continue 
90 KK = KK + J + 1 

Both segments of the computer code are of interest, The first part, with 
= 0 can be applied to triangularize a system of Npx + m equations in Npx 
variables; the latter part can be used, with Npx = 0, for recursive triangulariz- 
ation of an overdetermined system where the results are stored in vector form 
to economize on storage. 


II. FILTER: TIME UPDATING (FROM T TO T + DT) 

The time updating uses an efficient and compact algorithm that utilizes p 
only one component at a time; i. e. , by considering n^ separate mappings 


X k+1 * X k + V k)P " <k) ) 

| k = 1, • • • , n p (B-l) 

p n+l (k) = M k p n <k) +w n <k) ) 


where 


X 1 = X n andX n +1 = x n+l ’ see Ec *- 3 "4. 
P 
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The data equation for smoothing corresponding to Eq. 3-12 turns out to be 


£ 

j=l 


v (k)p (i (k) ^ s’* (k . j) p n (k + j) 


+ ^ s' (k, 

i = n-K t 1 

» n 


j)p n + , ( k + j - n p ) 


Npx 

♦ £ " 


S (k, j) X k + ,(j - n ) 


j=Np+l 


y 

5 


♦ > R py (k ’ j> - *p‘ k » - •>» 


k s 1 , • 


n 

P 


(B-2) 


Inputs 

n p, n x , n^,, n d : are the dimensions of p, x, y and the number of dynam- 
ically occurring p components. 


T(n p ): Time constants for the colored noise parameters 
V <n x , n^): The first n d columns of the V p matrix correspond to the 
dynamic parameters. The last n p - n d columns are omitted because they are 
in theory zero. 


R (n ): Process noise standard deviation reciprocals, 

w p 

S(Npx + n p , Ntot): Ntot = Npx + + 1; the bottom n p 

used to store the smoothing related data equation terms. 


rows of S will be 
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Outputs 


S(Npx + rip, Ntot): The time updated S information array occupies the 

upper Npx rows of S, and tl,e lower n rows involve data equation information 
>!< P 

for smoothing, x (n ): Vector of smoothing related output terms, 

P 


DO 100 J = 1, Np @Np = n p 

If [t(J)*z] D(J) = EXP[-DT/r(J)] @z = zero 

100 Continue 

c Colored noise transition terms are set 


DO 106 JL = 1, Np 
If (Jl > Nd) Go to 102 
DO 101 1=1, Npx 
DO 101 K = 1, Nx 

101 S(I, 1) = S(I, 1) - S(I, Np +K)* V p {K, Jl) 

102 X = -R w (Jl)*D(Jl) 

x = X **2 

DO 103 K = 1, Npx 
w(K) = S(K, 1) 

103 x = x + w{K)s!(f|<2 

x ~ 'v/x 

X = X - x 

x' ,c (Jl) = x 
x = one/(xs;<\) 

DO 105 J2 = 2, Ntot 
6 = z 


@Nd = n^ 
@Nx = n x 


@w(Npx): work vector 


@one = 1. 
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c If (J2 = Ntot) 6 = \>!<55 w (Jl) 

c Colored noise is generally zero mean; thus z w is zero 
c in most cases. 

DO 104 I = 1, Npx 

104 6 s6+S(I, J2)*W(I) 

6 - 6y,«cr 

L = J2 - 1 

If (J2 > Np) L = J2 

S(Npx + Jl, L) = 6*X. ©Used for smoothing 

DO 105 I = 1, Npx 

105 S(I, L) = S(I, J2) + 6#w(I) 

c S(Npx+Jl, Ntot) = S(Npx + J1, Ntot) + &si<z^(Jl) ©Used for smoothing 
c 

6 = \>^R (J 1 )*<r 
w 

S(Npx + Jl, Np) = R (J 1 ) + 6 j;«X. @Used for smoothing 

w 

DO 106 I = 1, Npx 

106 S(I, Np) = 6*w(I) 


HI. SMOOTHING 




Suppose that smoothed values for the computed estimate, X * computed 
covariance, P c , and the sensitivity, SEN , corresponding to the augmented vector 
X = ^Jhave been obtained at time T + DT. Smoothed values for these terms, 
at time T, are obtained in the following way. Start with the previously stored 
smoothed array jV v (Np) j S >|! (Np, Npxy) [ z*(Np)J , (cf Eq. B-2), where 
Npxy = n + n + n , and the transition elements V (n x , n^). To facilitate 

p y . p 

algorithm implementation the columns of Eq. (B-2) are rearranged so that the 
terms p R (n referring to time T) displace those for P n ^! i» e., 
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(B-3) 





(k,n p ), s' (k, 1 ' r S^^k, n p 



k = 1, * " ,n. 


This arrangemont is bettor suited for the backward computation of 
P ( n p)' ' ’ *»P '(0* The backward smoother computations which are documented 
more thoroughly in Ref. 1 are summarized as follow^: 

For k = n , * * * , 1 recursively cycle through Eqs, B-4 through B-15: 

P 


Y: = 1 / «■' <k) 


{B-4) 


v(j); = s’^k, j)y , j - li * * • ,Npx 


(B-5) 


N 


px 


6* = v( j )X l' (j) 


(B-6) 


If (k < n d ) 


X*(k) = z^(k)y - 6 
P 


x c(" p + j) : = <(" p + j)- V j ' k,x * (k) 


(B-7) 


j = Nx (B-8) 


At the end of this sequence of n steps, X at time T + DT has been 
>;< P c 

replaced by X at time T. Compare this with Eq. 4-7. 


VyU): = Rpy(k.j)Y, 


j = t. 


(B-9) 
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6: 


N 


px 


» ^ v(j )SEN ''(j, m) + v y (m) 

J“1 


SEN* (k, m): = -6 


m = l f • * • , n 


SEN*(n +j,m): sSEN (n +j,m} 
P P 


+ V (j,k)6, j = 1 , * • ■ i n 


Compare this with Eq. 4-6, 


(B-10) 


N 


px 


s(j): 


E 

a= 1 


P*(j. •>*(«) . 

c 


j = 1, Npx 


(B-il) 






Of the equations involving P only this one involves the lower part of P , 

^ >|C C 

In the computer implementation only the upper part of P is used, 

c 


P ' (m, k): = z(m) 

v 


^(m.np+j); = P c (m ( n p +j) J j = 1, n. 


m = 1 , • * • , k - 1 


-a(m)V (j,k)\ (when k < n d ) 


(B- 12) 


N 


px 


P c (k,k): = v' 


y, v(j)zo) 
j=i 


(B- 13) 
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(B- 14) 


P> C V> j): = *s(j> j = k + 1, ’ ’ ' . n p 

P c( k ’ n p + a ) : = z ( n p + a ) " P c (k ’ k 'V ff ’ ^ 

p c( j>n P + ") ; = p c ( j ' n p +ff ) " z (i) v p ( ff > lc )» 

j s k+1, - • n p 

p c( n p + j’ n p +0, ) : " P c( n p + j' n p + ") 

- V j ' k)P c( k, ' l p + “) 

- z ("p + j) V p ( "' k) 

j = « 

(Of course, in the computer implementation, Eq. B-15 must be modified when 
k > n i, because V is only implicitly zero for these values of k. ) 

Cl p 

The FORTRAN mechanization which follows is, in fact, the mechanization 
currently employed in the JPL Orbit Determination Program. The apparent 
complication of the algorithm Eqs. (B-4) - (B - 1 6 ) belies the compactness and 
efficiency of the FORTRAN realization, 

Inputs 

X c (Npx): smoothed computed estimates of p and x (independent of y). 

P 1 (Npx, Npx): smoothed computed error covariance for the augmented 

c ,|{ _ >;« 
vector X . The lower triangular portion of P is not explicitly used in any of 
c c 

the computations. 
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SEN ' (Npx, Ny): smoothed sensitivity matrix (Ny = n^). 

Y (Nx,Nd): Transition matrix elements corresponding to the Nd dynami- 
cally occurring colored noise variables (Nx = n^, Nd “ n^). 

cr’(Np), S (Np,Npxyl): Data equation coefficients that were stored during 

the time propagation portion of the filter algorithm. Npxyt = Npx + Ny + 1 ; the 

»;< >!< 

last column of S is z . 

P 

X' P and SEN came from the previous smoothing step (or the filter 

c ® >|< ❖ ... 
values if this is the first step). V p , <r and S were stored during the time 

propagation portion of the filter algorithm and are now read back. 


Smoothing Outputs : P^, SEN * , X* 


DO 200 K = Np, 1, -1 

>:o;! >!t il"!! >;t :|c t;! j;t )!< »;! ;:t :l! 


102 


104 


106 


L = Np - K 

If (L = 0) Go to 110 
DO 102 J = 1, Np 
v(J) = S*(K,J) 

DO 104 J = 1, K 
S*(K. J) = v(J + L) 
DO 106 J = 1, L 
S*<K, J i K) = v(J) 


@ >;o!: 

(5) ^ Eq . B-3 


>1; s|; :;c s;: ;;o:t d! :lt ili * :;t II: S',! :l! 5i; sic :;t ;;i Jl: ;!; 


110 v = 1. /<r*(k) 
SUM = 0. 


@Eq. B-4 


DO 115 J = 1, Npx 

S*(K,J) = S*(K, J)s;; Y 


@Eq. B- 5 
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115 SUM = SUM + S*(K, J)>;«X*(J) 

X*(K) = S # ' (K, Npxyl )*Y - SUM 
If (K > Nd) Go to 125 
DO 120 J = 1, Nx 

120 X*(N +J) a X*(N +J)- V (J, K)>!(X " ! (K) 
c p c p p c 

125 If (Ny = 0) Go to 150 


@Eq, B-6 
@Eq. B-7 


DO 130 J = Npxl, Npxy 
130 S*(K,J) = S ' (K, J)>!<y 

If (K > Nd) Go to 150 


@Eq. B- 8 

@Npxl = Npx + 1 
@Eq, B-9 


c ;|t S|: * * ;Jt !jt * * s|; s|t * :|c * * ft* * :j: * i;s * * ❖ tj: :|t * 5|; :J; :I: :|t :!< * sjt >J: ;|t >;< :[( s|< * :|t * si; * * :|t 


135 


DO 140 L = l, Ny 
SUM = S *(K, Npx + L) 

DO 135 J = 1, Npx 

SUM = SUM + S*(K, J)*SEN*(J. L) 

SEN*(K,L) = -SUM 

DO 140 J = 1, Nx 


@**>:<* Eq. B-10 
(a)*** 


SEN*(Np + J, L) = SEN*(Np + J, L) + V (J,K)*SUM@* 

P 

.1. .i. .i. ,i. .i. .1. j. . i . .L .1 . .f, ,1, ,li *>, .t< kti •<< .)< ilj tli >l« tli »l» tli »►»».>■>»■» t*i jImN »lf 

i\i ;|s >ji sjs iji fj! :,c «]: #{« },» i,' JJ? >i» V s ( c *,> #,s * v vv v *i* »i» »»' 'i* v } i { *r •>* *r 


150 If (KS = FALSE) Go to 200 @KS = True (nominally) 

t'i 

c KS = TRUE if is to be computed 
DO 160 J = 1, Npx 
v(J) = 0. 

DO 155 L = 1, J 

155 v(J) = v(J) - P^(L, J)s;(S ’ (K, L) 

I = J + 1 

DO 160 L = I, Npx 
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v(J) = v(J) - P*(J, L)*S*(K,M 
c 

c 

>;< 

c 160 Loop in *£q. B-ll, using onLy the top part of P 

C 

c 


JJ = K - 1 
If (JJ = 0) Go to 180 
DO 170 L = 1, JJ 
P^L.K) = v(L) 

If (K > Nd) Go to 170 
DO 165 J = 1, Nx 

165 P*(L,Np + J) = P*(L,Np + J) - v(L)*V (J,K) 
c c p 

170 Continue 

180 P*(K,K) = y**2 

c 

DO 182 J = 1, Npx 

182 P*(K,K) = TP* (K, K) - S*(K. J)«v(J) 

c c 

If (K = Np) Go to 186 
JJ « K + 1 


(5)*** Eq. B-12 

@y,< 


@Eq. B- 13 


184 

c 


DO 184 J = JJ, Np 

P*(K,J) = v(J) @Eq. B- 14 


5|i >J< s|< >!■ >|s >J{ >|c >;< >]i sjs sjt :|c :J« :|c sj: !|t j|( >;< i|; :|s sjs ;J: i|s s!< s|( s|t >;« >[( :|i J|: ;Ji >!< >;« >|t >!:>!: :|s :[s 


186 DO 198 L « 1, Nx 
H = Np + L 
P*(K, U) = v(J) 

If (K > Nd) Go to 198 
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If (K = Np) Go to 194 
DO 192 J = JJ, Np 

192 P*(J, IJ) = P*(J,IJ) - v(J)*V (L, K) 

c c p 

194 P*(K,IJ) = P*(K,IJ) - P*{K,K)*V. (L,K) 

v* C C p 

DO 196 J = 1, L 

196 P*(N P + J,U) = P*(Np +J,IJ) - V (J,K)>;<P*(K,IJ) - v(Np + J)*V (L,K) 
c c pc p 

198 Continue 

q :[< >!t j[( >;< f|< ils i|! )|f iji j|: sj; j|t :|: ft :;e i|t >|t *’t >[; iji i|( tj: )|i :J; # # r;t :|s !|t ij! i|i >Jc s|i 

c The 198 loop is Eq. B-15 
200 Continue 

Smoothed estimates and covariances and smoothed consider covariances 
can be obtained by applying Eqs, 4-10 through 4-12. 
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